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I. INTRODUCTION 


Nonlinear transient analysis of structures, especially those involving 
material non linear! ties, are highly complicated even for the uniaxial case 
with a simple constitutive model. In the case of the beam bending 
problem for instance, there is a single component of strain 
along the longitudinal beam axis. This axial strain component is assumed 
to vary linearly over the cross-section by virtue of the assumption that plane 
sections remain plane. However, since the stresses are no longer proportional 
to strains in the inelastic range the distribution of stresses can be quite 
complex over the cross-section and along the length of the beam. Thus, for 
a beam which has yielded in a certain region the line of zero stress or strain 
does not coincide with the centroidal axis in that region. Hence, in general 
the centroidal axis is not an axis of constant strain for deformations in the 
inelastic range. 

A two-noded beam column element is often used to analyze inelastic 
response of models built up from such elements in conjunction with other 
elements. A linear axial displacement field and a cubic transverse dis- 
placement field are used in arriving at the stiffness properties of such 
an element. If the centroidal axis is used as the reference axis, it 
is not possible to satisfy equilibrium in the inelastic range. The 
same is true of the linear elastic range, if an axis other than the 
centroidal axis is used as the reference axis. Hence, a simple linear 
elastic analysis with an arbitrary reference axis can be used to demon- 
strate the point under consideration. 

Another feature of the nonlinear analysis by the finite element method 
is the necessity of integrating complex distributions of stresses and 
strain energy densities over the volume of the element. Because of the 



complexity of the integrand and the domain of integration, recourse to 
numerical schemes like Gaussian quadratures, Newton-Cotes, etc. has to be 
made in order to obtain approximate estimates of stress resultants and 
total energies. The number of integration points used over the cross- 
section and over the length of the beam element determines the degree of 
approximation. Accordingly approximate estimates of energies yield 
solutions of varying degrees of accuracy which depend upon the order of the 
integration scheme used. For an assessment of the quality of such 
approximate solutions it is necessary that an exact solution to a problem 
be known. Such a problem along with its accompanying exact solution has been 
outlined and a rigorous evaluation of the sensitivity of the quasi -static 
response to the order of the integration scheme is made. For problems, 
especially those involving transient response, for which no exact solutions 
can be obtained, only qualitative estimates of sensitivity can be obtained 
by a comparision of numerical solutions using different orders of integration 
schemes . 
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II. BEAM- COLUMN ELEMENT 


In this section the two-noded beam-column element of Fig. l.ai will be 
shown to yield inaccurate results when used to model a beam bending elasti- 
cally about its centroidal axis, but analyzed with respect to another 
reference axis. With the finite element formulation in mind the present 
discussion will be confined to Euler-Bernoulli beams subjected to concen- 
trated shears and moments. For such a beam the total potential energy of 
deformations is given by 
z 

= If I ^ - PiW(0)-P2w(<i)-Mt ^(0)-M2^U) 

0 A 

( 1 ) 

where u and w are the axial and transverse displacements of the reference 
axis and z is measured normal to the .reference axis. Minimization of tt 
with respect to u and w yields 

El ^ - EZA ^ = 0 (2) 

dx^ dx'^ 


EA ^ - EZA ^ = 0 C3) 

dx'^ dx^ 

as the governing equations (see Appendix A for details). Elimination of u 
from the above two equations yields the well known equilibrium equation 


El 


4 

d w ^ 0 


" dx^ 


(4) 


where 


I =I-7^A 
c 

in Eq. (5) is the moment of inertia of the cross-section about the 
centroidal axis which is separated from the reference axis by 2 . 
Equation (3) integrates to 

“ ° ^ ^ <=2 


(5) 


( 6 ) 
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Assume such a beam is to be modeled using the beam-column element of 
Fig. l.a. This finite element assumes a cubic transverse displacement 
of the reference axis. For satisfaction of equilibrium, Eq. (6) 
requires that the axial displacement vary quadratically over the length 
of the element. Hence, equilibrium cannot be satisfied by prescribing 
a linearly varying axial displacement field as is done in the case of 
the two noded beam element. As an example, a two noded beam element. 


used to represent a cantilever loaded at its free end with a concentrated 


transverse load P, yields the following results (see Appendix B for 

details): ^ 

U2 = 



(7-a). 

(7-b) 


^2 = - (3 + — — 

^ 1 OCT W ^ T 


dx^ X=£ 


1 ) 

c I 


(7-d) 


These results are correct only for the case I =1^ i.e., the reference 
axis for the two noded beam-column element is the centroidal axis. 

On the other hand the same cantilever if modeled with a single beam 
element of Fig. l.b yields the correct results and ensures a complete 
satisfaction of equilibrium in keeping with the elementary beam theory. 

From this linear elastic example it can be concluded that a three noded 
beam element is necessary to analyze the inelastic response of structures 
built up using frame elements, since in general the reference axis is not 


the neutral axis or an axis of constant axial strain when yielding takes 
place over portions of the element. 
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III. SENSITIVITY OF RESPONSE TO APPROXIMATIONS OF ENERGY ESTIMATES 


a. Exact Solution : 

The exact solution of a quasi-static inelastic beam problem will 
be developed in order to demonstrate the sensitivity of the response to the 
order of the integration scheme used in approximate solutions of the 
problem. Consider the cantilever beam of rectangular cross-section of 
Fig. 2. a subjected to a bending moment at its tip. The material of the 
beam is assumed to be linearly elastic — linearly strain hardening as 
shown in Fig. 2.b. For monotonic loading such a material is conservative 
and hence the principle of the stationary value of the total potential energy 
can be used to obtain the governing Euler-Lagrange equations for the beam. 

Euler-Bernoul 1 i hypotheses imply that 


(X) = -y 


dx‘ 


From Fig. 2.b 


a = E-, e 


= + E2(e - Ey) 


f- 


if e < e 

- y 


if e > £. 


The total potential energy of the beam of Fig. 2 is given by 


.-/w,dv ^/W2dv-H^l 
V, Vo 


where 


W, = 

I O 


1 '2 
2 


( 8 ) 

(9-a) 

(9-b) 

00 ) 

(11-a) 


5 


(n-b) 


W2 = ^ ^ + Oy(^-Ey) 


V^ is the volume of the beam within which the beam is everywhere elastic and 
^2 is the volume within which the strain at every point exceeds the yield 
strain, e^. 

From purely statical considerations it is obvious that the beam of 

d2 

Fig. 2 experiences a constant curvature, — ^ and, hence, d« of Fig. 3 is a pure 

constant and not a function of x. This can be shown to be true by minimizing 

the total potential energy with respect to»not only w,but also d 2 » assuming 

apriori that d 2 is a function of x. The Euler-Lagrange equation resulting 

from the variation of tt with respect to w and d« implies that the curvature, 
d^w 


dx 


2 ^»and d 2 are both constants. 

The two Euler-Lagrange equations 



( 12 ) 


^=0 (13) 

dx^ 

and the associated boundary conditions 

w(0)=^=^%i = 0 (14-a) 

dx'^ 

and 

(Eg-Ei )y (d2-d2)4|^[E,d|+E2(d^-d^)][^^4^>M (14-b) 


are the result of the variation of it with respect to d 2 and w. 
Equations (13) through (14) are satisfied by 
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(15) 


w=ax 


2 


where 


[(- 




E<i£ bd, 
1 y 1 


) + (1 . ^){l - (^)2j] 


[(J)^ + (^) {1 - (^)^] 

d] 


(16) 


Use of Eq. (15) into Eq, (12) yields 


adg = - ^ 


Cl 7) 


Substitution of a from Eq. (16) into Iq. (17) yields the cubic equation 


rp'^-3(r+s)p-2(l-r)=0 


(18) 


with 


p=A), and 


*1 


1 


El^ybd^^ 


(19) 


For any given values of r and s, Eq. (18) can be solved by trial and error 
to obtain the corresponding value of p. The total strain energy of deformation 
is then given by 

[{p^+(1-r)0-p^)l(2ad,)2-3reyCl-p^)C2ad,)-3e2p(^l.p)] (20) 

If the beam is unloaded after being loaded into the inelastic range and if 
it is assumed that unloading takes place elastically then it can be easily 
verified that the recoverable energy of deformation is given by 


•4^3 


2 . 2 2 


U^=E^bnd^[|p'^(ad,)‘+rSy‘(l-p)-2ad,eyrO-r)(lV)+ f(l-r)^(l V) (ad^)] (21) 
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and the energy dissipated in the process of unloading is 

U(,=Ut-Ur (22) 

b. Approximate Numerical Solutions : 

Three numerical solutions to the problem just outlined are obtained via 
the minimization of the total potential energy using the three noded beam 
elements of Fig. l.b. and using numerical integration schemes of three 
different orders for the computation of the strain energy of deformations [1]. 
The first two schemes use Gaussian quadratures of two different orders in 
that the beam is divided into two strips in one case and into four strips 
in the other case. A two point Gaussian quadrature formula is used in 
each direction (the length and the breadth of the strip). The third scheme 
is the Newton-Cotes scheme which uses a cubic interpolation with four 
points in each direction. The details of these schemes can be found in 
any book on numerical analysis [1]. It must be emphasized however that 
among integration schemes which imply that the integrand can be approximated 
by a polynomial of some order the Gaussian quadrature scheme is the most 
efficient. It is well known for instance that only n Gauss points are 
sufficient to integrate a polynomial of (2n-l) degree exactly while (n+1) 
points are necessary to integrate a polynomial of n-th degree exactly with 
the Newton-Cotes integration scheme. 

Table 1 shows a comparison of the exact response with that predicted 
fay the three different integration schemes. Of main interest is the dissipative 
energy. It is obvious from these results that a very high order quadrature 
is necessary for a good correlation with the exact solution especially for 
loads where yielding extends over a major portion of the cross-section of 
the beam. The higher the value of s in Table 1 the more extensive the 
yielding. 
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Such large variations in the estimates of deformations and energies for 
a quasi-static case suggest that a similar study be made in connection with 
the nonlinear transient response of structures. Since exact solutions to 
the nonlinear transient response of structures are virtually nonexistent one 
has to rely on a comparison between numerical solutions using different 
integration schemes. 

Figure 4 shows the details of a solid cross-section beam clamped at both 
ends and loaded impulsively at the center by a triangularly varying impulse. 

The material characteristics of this beam are assumed to be as shown in Fig. 2 
with 0 ^ = 28.96xlO^N/m^ and = 68.95xlO^N/m^ and (E 2 /E^ )-0. 4286. A numerical 
analysis of the transient response of this beam is performed using direct energy 
minimization at each time step. For the purposes of the integration of strain 
energy the two different Gaussian quadrature schemes shown in Fig. 4d are employed. 

Figures 5 and 6 bring out quite vividly the differences in the response 
resulting from the use of the two integration schemes. It must be remarked 
that the deformation- time plots for the two schemes are very nearly identical 
and hence are not shown. The same cannot be said of the acceleration or 
energy plots, however. As regards the acceleration- time plot, at certain 
instances the accelerations predicted by the two models can be seen to be off 
by as much as 100%. The nonlinear character of the model coupled with num- 
erical approximations defies predictions of trends in response as evidence 
by the crossing of the two response curves. The model with 16 Gauss points 
can, in general, be expected to provide a better estimate of dissipation 
than the model with 8 Gauss points. It can be seen from Fig. 6 that the model 
with 16 Gauss points predicts a higher dissipation than the model with 8 
Gauss points by as much as 15% over the interval considered. Again, this 
seems to be a peculiarity of these two models and in general one model which 
may at any given time predict a dissipation higher than the other may very 
likely also predict a dissipation lower than the other at another time. 
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Next the sensitivity of the response of thin-walled frame elements is 
examined. Such elements can be expected to "tone down" the sensitivity of 
response, if any, to the order of the integration scheme. Figure 7 shows half 
the finite element model of a horizontal box beam supported on four uprights 
and loaded impulsively at the center. The beam and the supports are assumed to 
be made from the same material the stress-strain curve for which is similar to 
that of Fig. 2. For the purposes of numerical integration using Gaussian 
quadratures each wall of the box section is assumed to be divided into 
rectangular strips extending between the two nodes of a frame element with 
six degrees of freedom at each node. Table 2 lists the relevant properties 
of the model shown in Fig. 7. Two different integration schemes are employed - 
one in which each wall is divided into two strips and another in which the 
same is divided into four strips. As before a two point Gaussian quadrature 
formula is used in each direction. In the analysis, the effects of shear 
deformations are ignored in the interest of simplicity of the constitutive 
model in the inelastic range. An impulsive load is applied at node 3. 

Figure 8 shows the plots of acceleration versus time and Fig. 9 shows the 
plots of dissipative energy versus time for the two integration schemes. 

Although the response is very nearly identical in the initial stages the 
responses for the two cases diverge from each other significantly with the 
passage of time. As expected the structure is much less sensitive to the 
order of the integration scheme than the solid beam of Fig. 4. However, it 
is again evident that no trends can be determined since both integration schemes 
are only approximate. In fact, it would seem that as a result of the approxi- 
mations in calculations of the strain energies of deformations upper bound 
solutions, even for conforming finite element models, are not guaranteed. 

Thus one scheme may provide a better answer than the other at any given time 
but the higher order scheme may be expected to give better results overall. 
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IV. CONCLUSIONS 


This study has established the necessity of examining every finite 
element used for modeling structural behavior in the inelastic range, for 
consistency of assumptions that v/ill guarantee satisfaction of equilibrium. 
Furthermore, this study has revealed that structures consisting of frame 
members wherein numerical integration has to be used for the purposes of 
evaluating stress resultants or strain energies will be sensitive as regards 
their accelerations and energy dissipations to the order of the numerical 
integration scheme used. This sensitivity can be expected to be more 
pronounced for frame members v/ith solid cross-section in comoarison vnth 
thin-walled members. This study demonstrates the need for higher order 
integration schemes for improved quality of response in the nonlinear range. 
Such higher order integration schemes however, would only make the already 
expensive numerical analysis of nonlinear response only more so. 
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Table 1 : Comparison of Beam Response for Various Integration Schemes 


Scheme 


Recoverable Energy 

1 

Dissipative Energy 

s 


L 

^2 

^'i 

^2 

^1 

^2 

Exact 

1.122489 

1 

1.587835 

171.19222 

269.8933 

16.0936 

105.2099 

Gauss with 
32 stress i 

ref. pts/ 
el ement 

1.13313 

1 . 58667 

171.1086 

269.8933 

19.5555 

85.24097 

Gauss v/ith 
16 stress 
ref. pts/ 
element 

1.06848 

1.63969 

1 170.6667 

! 

267.6774 

0.5803 

109.20314 

Newton- cotes 
v/ith 32 
stress ref. 
pts. /el ement 

1.15733 

' 

1.53067 

1 

159.5500 

1255.6100 

37.0460 

87.7340 


v/here s ^ = - 0.8466, ~ " 1-0582 



































Table 2. Properties of the Beam on Elastic Supports 


Nodes 

m 

B 

Lumped Mass 
Ki logram 

Element 

Length 

m 

cm 

cm 

B 

1 

0.8001 

0.0 


1 

0.7620 

3.81 

1.463 

.0704 

B 

2.7051 

0.0 

- 

2 

0.7620 

3.81 

1.463 

.0704 

H 

0.0 

0.762 

0.07793t 

3 

0.40005 

2.9261 

7.62 

.141 

B 

0.40005 

0.762 

0.07793 

4 

0.40005 

2.9261 

7-62 

.141 

B 

1.7526 

0.762 

0.1852 

5 

0.9525 

2.9261 

7.62 

.141 


0.8001 

0.762 

0.14988 

6 

0,9525 

2.9261 

7.62 

.141 

fl 

2.7051 

0.762 

0.11138 







tA non-structural mass of 0.07793 kilogram is assumed to exist at node 3. 


H 















U(0=(1**?)U^+? U 2 , 5=X/A, 

w(c)=0-3c^+2c^)w^+£(C-2?^+5^)e^+(35^-25^)W2+5'("C^+e^)e2 
Figure l.a. Two Noded Beam-Column Element 



U(c)=(1-3C+25^)U,+(45-45^)U2+(-5+25^)U3 
W(5)=(l-35^+2s^)W^+t(5-25^+«^)e^+(35^-2c^)W2+«(-5^+5^)02 
Figure l.b. Three Noded Beam-Cttlumn Element 
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I = .508 m , b = .0254 m , 

^ £ = density = 0.00774 kg/cm 


a. A typicol Beam clamped at both ends 


b. Beam cross-section 


c. 



2.2241 X 10 Newtons 


t = 0 t = 5 X 10^ t= 1.0 xlO®(seconds) 


Variation of load P 


Izzd — ' Points 

Scheme I Scheme H 

d. Strain energy integration schemes 



Figure 4. A Solid Cross-Section Beam Under an Impulsive Load 





Acceleration x 10 m/sec 



Dissipative Energy x 10" Joules 







e 8. Accel 


'’.0 


21.0 


25.0 


29.0 


33.0 


Time x 10 


5 


Seconds 


ration-Time Signature of Node 3 (Fig. 7). 


Dissipative Energy x 10" Joules 



Figure 9. Dissipative 



J 1 1 i I I I 

25.0 29.0 33.0 

10® Seconds 


Versus Time (Fig. 7). 


APPENDIX A 


Consider an Euler-Bernoulli beam subjected to end shears and moments. 

For analyzing the response of this beam, its deformations are referenced 
with respect to an axis which does not coincide with its centroidal axis 
around which it bends. This axis is designated as the reference axis. 

The potential energy expression for such a beam is given by 

IT = / f - z A]^dAdx-P,w(0)-P2w(i)-M,^(0)-M2 (A-1) 

which upon performing the integration with respect to the area of cross- 
section reduces to 


IT 


EA 

2 


I 

f 

0 


(f)^dx-EZA 


/“(^)(^)dx + /(^)^dx-P,w(0)-P2w(O 

0 dx'^ ^ 0 dx'^ ‘ 


M clw(O) ^ » _^(a) 
'^l dx ^2 dx 


where ZA and I are respectively the first and the second moments of the 
area about the reference axis. 

Upon requiring that the variation, 6it, of the potential energy vanish 
the necessary Euler-Lagrange equations or the equations of equilibrium of 
the beam are obtained. Thus, 

aw=EA /(|j)a(^)dx-ElA /(^)5(^)dx-EZA /^(^)(^)dx 
0 ax OX Q ax q ax 

£-2 2 

+EI A (^)6(^)dx-P,Sw(0)-P26w(O-M^6{^(0)) - H2a(^(^))=0 (A-3) 
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After integrating by parts the following equations and boundary conditions 
are obtained. 


El - EZA ^ = 0 
dx^ dx'^ 


EA - EZA ^ = 0 
dx*^ dx'^ 


Either 


EA- 


du(0) _ d w(0)_ ^ 


dx 


dx 


- eza = p 

dx^ dx^ ' 


ei^M . eza ^ = 


dx 


dx 


- M. 


EA^^ - EZA 

dx'^ 


= 0 


- EZA ^ 

dx^ dx 


= - P, 


rr d^w(^) _ rjA du ( g. ) , 

. 2 dx 

dx 


= M, 


(A-4a) 

(A-4b) 

or 

u(0) = ug 

w(0) = wg 

^(0)_ * 
dx " 

(A-5a,f) 

u(5,) = u* 
w(!i) = w* 


The starred quantities in the above equations denote prescribed quantities. 
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APPENDIX B 


a. Cantilever Beam Modeled with a Two-Noded Beam Element : 

The two noded beam-column element of Fig. l.a. is used to model 
a cantilever beam subjected to a load P at its free end. The potential 
energy for such a beam is given by Eq. (A-1) with P^=M^ =1^2=0 and p 2 =-P. 

It can be easily verified that the substitution of expressions for u(c) 
and w(?) from Fig. l.a. into this potential energy expression yields 

+ C12w^+12w-j6-j£-24w-^W2+12w-| 02i?'+40^ji^-12w26iii 

2 2 2 2 

+40^02^1 + 12w2 -1 2w20 2^+402^2. ] + PW 2 (B-1) 

Upon requiring that the variation of tt with respect to the nodal variables 
^1’ ^2’ '^l’ ^^2* ®1* ®2 following equations are obtained. 

Either or 

f^(u,-U2)+ ^ ( 62 -e ,)=0 

f (u2-U,)+-^(6^-82)=0 U2=u| 

■^[12(Wi-W2)+6i,(e^-02)]=O w^=w^ 
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Either 


or 


^[12(w2"W^ )-6a(0-j+02)]+P-O 


Ui =Uf'^ 

”2 ”2 


™ (u2-u^)+^[6(w^-W2)+24(2e^+62)]=0 


(B-2a,f) 


e^=e» 


™ (u.|-U2)+^[6(w^-W2)+2ji(e,+262)]=0 


62=e| 


Since, for the cantilever beam the nodal displacements and 0^ 

are all prescribed to be zero the three equations corresponding to unknown 
displacements U2» W2 and ©2 simplify to 


U2=Z02 


El 


^(12w2-602^-)+P=O 


(B-3a,c) 


U2+^(-6W2+402O=O 


The above three equations when solved simultaneously yield the following 
results 

2 


Uo= 


.. - P^~ C\ 

«2 “ 121 


I 


(B“4a,c) 


0 = - ^ 
®2 2EI 
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where 


Ic=i-az2. 
d^w 


Furthermore, — 5- at the free end (i.e. at x-a) which is proportional to the 
dx^ 

bending moment at the free end can now be evaluated to be 


2 T 

d vija) _ P5; (_c T \ 

^2 " ■ 2EI M ’ 


(B-4d) 


Clearly, the above results reduce to the correct strength of materials results 
if 1^=1 which is to say the reference axis coincides with the centroidal axis. 


b. Cantilever Beam Modeled with a Three-Noded Beam Element : 

The expression for the potential energy in this case becomes 


" = It W3- F ¥- 


EZA 

2“ [4u^w.|“8u2Wi+4u2W-j+3u.je-j2-4u2eiJi+U26-jJl-4u^W2+8u2W2 


- 4 U 3 W 2 +U 1 e2^-4u202^‘'‘3U3e2^] 

+ [12w^+12w.j0^£-24w^W2+12w.|02Ji+40^Jl^-12w2e-j£+40.j02J^^ 

+12W2-12W202S' +402^^] + PW 2 (B-5) 

Upon requiring that the variation of tt with respect to the nodal variables 
u^, U2S u^* w-j, W2 and 0-j, 62 be zero the equations corresponding to the 
unknown variables U2» u^* ^2 
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^2* ~ ■^^(8W2-402^)“O 


8 . X EZA, 


^3* I ‘^2'^ *^3^“ ■^(-4w2+302^t)-O 


EZA. 


^2* " ™(8u2-“4U3)+^(12w2-602il)+P-O 


(B-6a,d) 


02: - {-4u2+3u3)+^(-6w2+402S-)=O 


El 


Simultaneous solution of these equations yields 


^3 = 




2EI 


^2 " 


-z (M-) 

c 


(B-7a,d) 


W2 = 


Pr 

3EI, 


0 = - 
®2 2EI 


These are identical with the well known strength of materials results. 
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